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Abstract. We present a novel, fast and precise method for including the effect of light 
neutrinos in cosmological iV-body simulations. The effect of the neutrino component 
is included by using the linear theory neutrino perturbations in the calculation of the 
gravitational potential in the A^-body simulation. By comparing this new method 
with the full non-linear evolution first presented in [1], where the neutrino component 
was treated as particles, we find that the new method calculates the matter power 
spectrum with an accuracy better than 1% for "^rrii, < 0.5 eV at z = 0. This error 
scales approximately as (^m^)^, making the new linear neutrino method extremely 
accurate for a total neutrino mass in the range 0.05 — 0.3 eV. At z = f the error is below 
0.3% for '^rrii, < 0.5 eV and becomes negligible at higher redshifts. This new method 
is computationally much more efficient than representing the neutrino component by 
A''-body particles. 
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1. Introduction 

Next to photons neutrinos are the most abundant particles in our Universe. The 
fact that at least two of the three neutrino mass eigenstates have masses much 
larger than the current temperature, rrii ^ Tq, means that neutrinos contribute to 
the matter density and are important for cosmological structure formation. The 
mass differences established by oscillation experiments, Am^g — ^ 10"^ eV^ and 
|Am|3| ~ 2.4 X 10"'^ eV^ (see e.g. [2] for a recent data analysis), mean that the heaviest 
mass eigenstate must have a mass of at least m ~ 0.05 eV. This in turn will have the 
effect of suppressing the power spectrum of matter fluctuations by ~ 5%, an effect which 
is significantly larger than the precision with which the matter power spectrum can be 
measured in upcoming surveys. 

This has two interesting consequences, first it means that neutrino mass must 
be included in the analysis of future data in order to avoid seriously biasing the 
measurements of parameters such as the dark energy equation of state, w [3]. Second, 
it means that future large-scale structure surveys potentially have the power to probe 
neutrino masses as low as the current lower bound from oscillation experiments [1] . 

Consequently this has led to a much increased interest in gaining a detailed 
understanding of how neutrinos affect structure formation. In linear theory the effect is 
extremely well understood and can be calculated to a precision much better than 1%. 
However, these results apply only on very large scales, k -C 0.1 /iMpc~^, whereas most 
of the cosmologically relevant information from large-scale structure surveys is in the 
range k ~ 0.1 — 0.5/zMpc~^. On these intermediate scales it is mandatory to correct 
for non-linear effects, even if surveys are carried out at intermediate redshifts where 
non-linearity is weaker. 

Non-linear corrections can be found in a number of ways: The most accurate, but 
also most time-consuming method is to use full A^-body simulations with neutrinos 
included in a self-consistent way. In a previous publication we carried out a large suite 
of simulations with neutrinos included [1] , and found a significant non- linear correction 
caused by neutrinos. The effect is large enough to be important even for masses as 
low as 0.05 - 0.1 eV. Another way is to use higher order perturbation theory which, 
however, is applicable only up to /c ~ 0.1 — 0.2 /iMpc^^ today, if the precision needs to 
be better than 1-2%. Finally it is possible to use semi-analytic methods such as the halo 
model which has been calibrated using ACDM simulations. Neutrinos can be included 
fairly easily in this scheme, but the accuracy can only be tested against full A^-body 
simulations. 

In conclusion, a very large number of high-resolution A^-body simulations with 
neutrinos included will be necessary for future data analysis. The method used in pQ 
is well suited for large neutrino masses with ^mi, > 0.5 eV at low redshift, but for 
all masses the simulations inevitably become noise dominated at high redshift, unless 
extremely many neutrino A^-body particles are used. 

In the present paper we describe a novel, fast and precise method for simulating low 
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mass neutrinos in A^-body simulations which uses the fact that low mass neutrinos have 
very little non-linear clustering even at low redshift. Instead of simulating neutrinos as 
particles with individual velocities we describe the local neutrino density on a grid. This 
density is evolved forward in time by using linear theory. The CDM/baryon component 
is followed simultaneously using the full TreePM method, with the neutrino component 
contributing to the long range force calculated using the PM method. As will be shown 
below the accuracy of this method is better than 1% for the matter power spectrum 
today on all scales for all neutrino masses below ^mj, ~ 0.5 eV. 

The method is very powerful, both because it speeds up simulations with neutrinos 
very significantly (by factors of order ~ 10), and because it can be applied to any other 
component of energy density which is almost linear. This could for example be dark 
energy with linear perturbations. The only requirement is that the linear component 
can be included in a Boltzmann solver like CAMB [5]. 

We note that our new method is conceptually similar to the perturbative approaches 
presented in [6],[7j, except for the fact that we give the CDM component a full non-linear 
treatment. 

In Section [2] we describe our cosmological model and the set-up of initial conditions 
for our neutrino A^body simulations. In Section [3] we present our results and in Section 
m we show that our results have converged. Finally, Section O contains a discussion and 
our conclusions. 

2. Initial Conditions 

2.1. The Cosmological Model and Particle Initial Conditions 

As long as the evolved primordial density perturbations set down by inflation remain 
small, their evolution can be calculated precisely using the linearised Einstein and 
Boltzmann equations [8]. However, once structure enters the non-linear regime precision 
studies require the use of A'-body simulations. To set up the initial conditions (ICs) for 
our simulations we have calculated the transfer functions (TFs) using CAMB [H]. The 
further evolution is followed in gadget-2 [9] run in the hybrid TreePM mode. 

We have assumed a fiat cosmological model with density parameters VL\^ = 0.05, 
Qni = 0.30 and Q\ = 0.70 for the baryon, total matter and cosmological constant 
components, respectively, and a Hubble parameter of h = 0.70. We vary the CDM 
and neutrino density parameters (f^coM and Q^,, respectively) such that they fullfill the 
condition ^cdm + = 0.25. We have assumed a primordial power spectrum of the 
standard scale-invariant Harrison-Zel'dovich form. The amplitude gives erg = 0.878 for 
a pure ACDM model. 

To generate the ICs for our simulations we have built a parallelized IC generator 
directly into GADGET-2. The TFs are used to generate the position and velocity ICs for 
the A^-body particles. In addition to the Zel'dovich Approximation (ZA) [10], we have 
included a correction term from second-order Lagrangian perturbation theory (2LPT) 
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[m [I2]. Due to our high A^-body starting redshift [z = 49), this second-order term is 
small, and therefore basically does not affect our results. 

The neutrino component is included in the A^-body simulation either as real space 
particles or on a grid in Fourier space. The CDM component is always treated as 
particles and for its initial power spectrum we have used a weighted sum of the CDM 
and baryon TFs. Gas physics is not included in the iV-body simulation, since it does 
not significantly affect the scales simulated. The initial conditions for the CDM and 
neutrino components have been generated with the same set of random numbers, to 
enforce the assumption of adiabatic initial fluctuations. In the neutrino particle method 
we add a thermal velocity drawn from a relativistic Fermi-Dirac distribution to the 
neutrino A^-body particles. When considering the power spectrum this thermal velocity 
acts as a wavenumber dependent suppression term (see [1] for further details). 

2.2. Initial Conditions and Evolution of the Neutrino Grid 

When the linear neutrino component is represented on a grid in the A^-body simulation, 
there are no neutrino A^-body particles. The representation is done as follows. A 
realisation of the neutrino TF is generated with the ZA only, using the same set of 
random numbers as was used to make the initial conditions for the CDM particles. The 
neutrino component stays in the Fourier domain throughout the A^-body simulation. 
The linear neutrino modes are added directly to their CDM counterparts whenever 
the long range force is calculated, taking proper care of the fact that only the latter 
modes should be deconvolved from the smoothing effect of the Cloud-in-Cell (CIC) mass 
assignment used in gadget-2. Whenever the long range force is calculated the neutrino 
component is evolved by interpolating between a library of TFs which span the redshift 
range over which the power spectrum is to be evolved in the A^-body simulation. 

The neutrinos are not included in the short range Tree part. To compensate for 
this fact the neutrino modes are not smoothed in Fourier space with the Gaussian factor 
exp(— /c^r^), where rg = 1.25 h~^Mpc, as are the CDM modes (see [13] for an analysis of 
the TreePM method). This lack of smoothing could be a problem since the accuracy of 
the PM method breaks down as the PM grid mesh size is approached (see Section Hj). 

A concern about the neutrino grid method could be that there are rotations in the 
CDM density field in real space, which could bring the CDM and neutrino components 
out of phase. But as long as only the linear neutrino density distribution contributes 
to the matter power spectrum these offsets due to rotations are negligible. Another 
more serious problem could arise since the linear part of the neutrino power spectrum 
and the CDM power spectrum are evolved in two different time integrators, CAMB and 
GADGET-2, respectively. Force calculation errors in the two integrators could lead to 
discrepancies, which would not only affect the absolute power spectrum but also the 
relative power spectrum between the neutrino grid and particle approaches. In Section 
Hlwe address these issues. 
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Table 1. Parameters for the iV-body simulations used to make the power spectra 
presented in this paper. A^cdm and iVi/,part is the number of CDM and neutrino TV- 
body particles, respectively. A^i/^grid is the size of the neutrino Fourier grid. ^ is 
the total neutrino mass, and it is in all cases related to the one-particle neutrino mass, 
TOi/, by = 3m,. r2,_o is the fraction of the critical density contributed by the 

neutrinos today. A^pm is the number of one-dimensional PM grid points and -Rbox 
is the size of the simulation volume. All simulations have a starting redshift of 49. 
Note that in model Cio the neutrino distribution has been kept totally homogeneous 
in the A^-body simulation. We will use the notation Xi/Xj to indicate a relative power 
spectrum calculated between the models Xi and Xj from (Xi/Xj — 1) ■ 100%. 

3. Results 

This section will focus on our main results, whereas the next section will contain a 
description of detailed convergence tests. Table [T] shows most of our performed A^-body 
simulations. 

Fig. [Ushows the difference in the total matter (neutrino plus CDM) power spectrum 
for 3 different neutrino masses, between simulations where the neutrino component is 
represented either on a grid or as A^-body particles. The difference between the two 
methods are shown at various redshifts. At our initial A^-body starting redshift, z = 49, 
the two methods are identical, as they should be. In the simulation where the neutrino 
component is represented by particles the neutrino A^-body particles free-stream out 
of the gravitational potential wells and generate a white noise term as soon as the 
simulation is started. This white noise term affects the difference power spectra in 
Fig. [T] for k > 0.1 — l/iMpc"^ (the actual wavenumber depending on redshift, neutrino 
mass and the number density of neutrino A^-body particles), and is solely due to the 
finite number of neutrino A^-body particles. That this is indeed a white noise term can 
be seen from Fig. [2] where the neutrino power spectrum is shown for all neutrino masses 
simulated. This white noise term can furthermore be seen in Fig. El where the neutrino 
A^-body particle density distribution is shown at z = 4 (i.e. the central panel). 

From Fig. [2] it can be seen that the neutrino A^-body particle white noise term is 
frozen on small scales, i.e. there is a maximum white noise level, which is determined 
by the number of neutrino A^-body particles. Note that for = 0.3 eV and z = 24 

even the neutrino fundamental mode is white noise dominated. This noise term persists 
and can be identified in the top panel of Fig. [H 

As the simulation evolves the CDM density perturbations grow so that the neutrino 
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Figure 1. Percentage differences in the total matter power spectrum at different 
redshifts between simulations where the neutrino component is represented in the N- 
body simulation either on a grid or as particles. A negative power difference indicates 
more power in the particle simulations. From top to bottom the total neutrino 
mass is 0.3 eV {B1/B2), 0.6 eV (C1/C3) and 1.2 eV {D1/D3), respectively. All the 
power spectra presented in this paper have been calculated on a 1024'^ grid, using a 
deconvolved CIC mass assignment scheme for the A^-body particles. 
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Figure 2. Neutrino power spectra for a total neutrino mass of 0.3 eV (left), 0.6 eV 
(middle) and 1.2 eV (right) at various redshifts. The linear theory neutrino power 
spectrum, convolved with our chosen random numbers, are shown with solid lines 
(models Ci and i?i), simulations with 256^ neutrino particles with dotted lines 
(models B2, C2 and D2) and finally simulations with 512'^ neutrino particles are shown 
with dashed lines (models C3 and D^). 

white noise term in the matter power spectrum dominates only on ever smaller scales. 
This is clearly seen in Fig. [T] where the agreement between the grid and particle 
representations improves on small scales as the simulation evolves. At 2; = 4 the two 
methods give basically identical results until k > l/iMpc"^. The disagreement beyond 
this wavenumber is caused by the white noise term. 

In sum, for '^m,y = 0.6 eV the neutrino grid method is accurate at the 0.1% level 
for all scales until z = A. This very low discrepancy is impressive given the fact that the 
neutrino component contribute up to more than 10% at this redshift (see the left panel 
of Fig. [H]), as well as the fact that the neutrino component is evolved in two different 
integrators. 

As the redshift falls below 4 the two representations begin to differ in the range 
k ~ 0.1 — 1 /iMpc"^. This difference is a non-linear correction coming from the fact 
that the neutrino A^-body particles are coupled to the non-linear gravitational potential 
whereas the neutrino grid is only evolved in linear theory. The wavenumber range where 
this difference appears can be explained by the convolution of two terms. A non-linear 
term growing rapidly on small scales and the fact that the neutrinos contribute most on 
large scales, as can be seen in Fig. [HI 

The lowest wavenumber at which this non-linear correction term becomes 
important, can also be identified in the middle and right panel of Fig. [21 as the range 
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Figure 3. Density grids for tlic CDM (left), neutrino particle (middle) and neutrino 
grid (right) components. In all cases the total neutrino mass is 0.6 eV. The top 
row is at z = 49, the middle row at z = 4, and the bottom row at z = 0. In the 
bottom row the square root has been taken of the first two density distributions. The 
images are centered at the highest density region in the simulation volume and they 
have a thickness of 20 /i^^Mpc and a side length of 512 ft,~^Mpc. The particle density 
distributions are found using the adaptive smoothing length kernel from [14] (taken 
from model C3), and the neutrino grid density distribution is an inverse FFT of the 
linear neutrino Fourier grid imbedded in the A^-body simulation volume (model Ci). 

from where the particle neutrino power spectra break away from the hnear theory 
evolution. The highest wavenumber at which the non-linear correction term matters 
is not only determined by non-linear neutrino modes at that scale but also by mode- 
coupling between the CDM perturbations at that scale and the extra non-linear neutrino 
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Figure 4. Evolution of tlic difference in the total matter power spectrum between a 
pure ACDM model (Ai) and a model with = 0.6 eV neutrinos on a grid (Ci). 

The difference expected from linear theory is also shown (black solid lines). 

contribution at larger scales. 

In the 'Y^mjj = 0.6 eV case at z = the two methods differ by at most 1.25% at 
roughly k ~ 0.25 /iMpc"^. Focusing on the = 1.2 eV cas^, we see from Fig. [1] 

that the non-linear correction is at the 5% level. As expected the non-linear correction is 
greater as the neutrino mass is increased because Qi, increases and the neutrino thermal 
velocity decreases. Since Qi, is proportional to the total neutrino mass and the thermal 
velocity is roughly inverse proportional to the one-particle neutrino mass this explains 
the {Y2m,y)'^ dependence on the size of the non-linear neutrino correction term. This 
scaling can also be seen to roughly hold in the ^ m,^ = 0.3 eV case, where the maximum 
non-linear neutrino correction to the matter power spectrum is at a negligible 0.25% 
level today (the discrepancy at large scales is due to a finite number of neutrino A^- 
body particles, see Section Hj). Note that as expected the wavenumber corresponding to 
the maximum non-linear neutrino correction propagates slighly from k ~ 0.3 /iMpc"^ 
in the ^rn^ = 1.2 eV case to /c ~ 0.2/iMpc~^ in the '^m^ = 0.3 eV case. In sum, 
the percentage maximum non-linear correction to the neutrino component today is well 

if This high neutrino mass was only used to illustrate when the linear neutrino approach breaks down. 
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fitted by the relation m,y/0.54eV)^. 

Finally, in Fig. H] we show the evolution of the difference in the matter power 
spectrum between a pure ACDM simulation {Ai) and model Ci with 0.6 eV neutrinos 
represented on a grid. For the scales simulated it can be seen that the turn-over in the 
difference power spectrum is created at low redshift, z ~ 1 — 4, and that it propagates to 
larger scales. The almost perfect agreement between linear and non-linear theory seen 
in the figure at z = 24 can in practice only be achieved by representing the neutrino 
component on a grid. 

4. Convergence Tests 

In the two approaches the neutrino component is evolved in different integrators. In 
the grid approach the neutrino component is evolved in CAMB with the full general 
relativistic equations and no cosmic variance. In the particle approach the neutrino 
component is evolved in gadget-2, with an accuracy limited by the finite box size, a 
finite number of CDM and neutrino A^-body particles (particle shot noise), as well as 
finite time steps and force resolution. It is important that the neutrino component is 
evolved accurately in both integrator^, so that our results are not compromised. 

Since we are interested in quantifying the non-linear neutrino correction term, we 
need to verify that the relative differences in the power spectrum between the two 
approaches have converged sufficiently. This demand is not as challenging as simulating 
a converged absolute power spectrum. We note that small-scale differences in the 
total matter power spectrum are acceptable since this part is dominated by the CDM 
component, which is evolved in the same integrator in both the grid and particles 
approaches. Small inaccuracies at small scales will therefore cancel out when comparing 
power spectra which only differs at the few percent level. But since the two neutrino 
representations are evolved in different integrators, we need a converged absolute 
neutrino power spectrum at the scales where the neutrino component contributes to 
the total matter power spectrum. 

Those test runs described in this section which are not listed in Fig. [1] are evolved 
in a 512 /I'^Mpc box and all have 256^ CDM particles and a 512^ PM grid. 

4-1. Initial Velocities and their Evolution 

We have found the gravitationally induced initial A^-body velocities by taking the time 
difference between two displacement grids centered around our A^-body starting redshift 
by some small ^z. We have run 3 simulations varying Az in the range 0.5 — 2, and 
found that our choose of Az ~ 1 has converged at a few tens of a percent. 

In GADGET-2 particle velocities are redshifted, not the momentum. This Newtonian 
evolution is very well justified for CDM particles but could pose a problem for low mass 

§ The accuracy with which CAMB calculates the linear power spectrum is far better than 1% on all 
scales. 
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Figure 5. The figure sliows tlie effect of redsliifting eitlier tlie Fermi-Dirac thermai 
velocity or momentum from z = 49 until z = 0, for 3 different neutrino one-particle 
masses. 

neutrinos, since at high redshift their Fermi-Dirac distributed thermal velocities can 
approach half the speed of light. In Fig. [S] we have calculated the cumulative neutrino 
thermal velocity distribution today in two different ways for our 3 neutrino one-particle 
masses simulated. Either by redshifting the velocity or the momentum from our A^-body 
starting redshift of 49. From the figure it can be seen that for our highest simulated 
neutrino mass of 0.4 eV, it is very accurate to redshift the velocity. In the 0.2 eV case 
redshifting the velocity or momentum gives the same result except for a slight difference 
at the high end of the velocity distribution tail. This small difference will not affect 
our results. In the 0.1 eV case this difference is slightly larger, but since it is mainly 
the neutrino iV-body particles with a thermal velocity drawn from the low end of the 
distribution which contribute to the matter power spectrum, redshifing the velocity is 
accurate enough. Likewise relativistic velocity addition is not necessary. Note that 
redshifting the velocity instead of the momentum leads to a lower thermal velocity 
and therefore more structure formation, so that from this point of view the size of the 
non-linear correction to the neutrino component in Fig. [T] is an upper bound. 

4-2. Box Size and Particle Shot Noise 

The finite size of the simulation volume as well as a limited number of A^-body particles 
can significantly affect the simulated power spectrum (see [15] for a recent analyses). We 
have chosen a box size of 512 /i~^Mpc from the considerations that a significant number 
density of neutrino iV-body particles was needed to suppress the thermal velocity white 
noise term on small scales. Increasing the box size by a factor of 2 would demand 
a factor of 8 more neutrino A^-body particles, stretching our computational resources. 
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Figure 6. Left: The non- linear neutrino correction as a function of the simulation 
volume and redshift. Dotted lines is for a 256/i^^Mpc box (C4/C5) and solid lines 
for our standard 512/i~^Mpc box (C1/C3). Right: Convergence as a function of the 
number of CDM A^-body particles at different redshifts. Dotted lines is for 256^ CDM 
particles (Ce/Cr) and solid lines for 512'^ CDM particles (C1/C2). In both figures 
Y^m^ = 0.6 eV. 




Figure 7. Convergence as a function of neutrino A^-body particles at various redshifts. 
Left: = 0.6 eV, dotted lines is for 256'^ neutrino particles (C1/C2) and solid lines 

for 512'^ neutrino particles (C1/C3). Right: '^m^ = 1.2 eV, dotted lines is for 256'^ 
neutrino particles {D1/D2) and solid lines for 512"^ neutrino particles [Di/D^). 



especially CPU time, to the limits. 

But, as shown in the left panel of Fig. [6l our chosen box size is sufficiently large for 
testing the neutrino grid method. The figure shows the non-linear correction from our 
highest resolution runs with ^m,^ = 0.6 eV (C1/C3), compared against a simulation 
where the box size, number of CDM and neutrino particles as well as the PM grid have 
been scaled down by a factor of 2 per dimension (C4/C5). Note that the two sets of 
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simulations have been started from different initial random numbers. The difference in 
the calculated non-linear correction is at the 0.2% level at most. Since we expect rapid 
convergence as the box size is increased, the chosen box size of 512 /i~^Mpc is sufficient. 

Including enough CDM particles is important for calculating the gravitational 
potential accurately since the neutrino iV-body particles move on this background. The 
difference between including 256^ {Cq/C-j) and 512^ {C1/C2) CDM particles on the 
calculation of the non-linear neutrino correction is seen in the right panel of Fig. [Ul At 
higher redshift the difference is largest at smaller scales, but the two sets of simulations 
begin to converge at low redshift. Note that the lower resolution CDM simulation 
predicts a slighly smaller non-linear neutrino correction at the 0.2% level. Again, 
expecting the size of the non-linear neutrino correction term to converge rapidly as the 
number of CDM particles are increased, our choose of 512'^ CDM particles is sufficient. 
Also note that the neutrino particle and grid methods approach each other at the largest 
scales as more CDM particles are included and that the differences between the two CDM 
resolutions mainly build up at high redshift, where the particle distributions are most 
homogeneous. 

Fig. [7| shows convergence as a function of 256"^ or 512"^ neutrino A^-body particles 
in the '^m^ = 0.6 eV and 1.2 eV cases. Focusing on the lower mass neutrinos it 
can be seen that the discrepancy between the 256^ and 512'^ simulations is largest 
at the smallest and largest scales simulated. At the smallest scales since here the 
neutrino white noise term contributes substantially to the total matter power spectrum, 
and at the largest scales since a noise term generated at high redshift is added to 
a real significant gravitational signal and therefore persists today. We conclude that 
in the range A; ~ 0.1 - 1 /iMpc"^ where non-linear correction terms to the neutrino 
component are important the difference between the neutrino grid and particle methods 
has converged as a function of the number of neutrino A^-body particles. 

In the = 1.2 eV case the noise term is even more pronounced at small scales, 

k > 0.3 /iMpc~^, since now ^l^ is larger even though the thermal velocity is smaller. As 
demonstrated in [1] the noise term decreases rapidly as the number of neutrino A^-body 
particles is increased, therefore we expect the difference between the neutrino grid and 
particle methods to have almost converged for 512^ neutrino A^-body particles. 

Since the maximum non-linear neutrino correction in the 1.2 eV case has almost 
converged for 512^ neutrino A^-body particles, and it has converged for 0.6 eV neutrinos, 
scaling to the 0.3 eV case we see that the 0.25% correction here should also have 
converged. 

With respect to the neutrino grid we (usually) have a Nyquist frequency of 
/iMpc"^ This is an oversampling of the neutrino perturbations, so that linear neutrino 
modes are not missing at scales where they contribute. We chose a cut-off in the neutrino 
grid modes at the grid Nyquist frequency. 
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4-3. Time and Force Resolution 

To get the neutrino Fourier modes at a specified redshift we have interpolated hnearly 
in the scale factor between the TFs in the library. These TPs are themselves equally 
spaced in the scale factor. To test the accuracy of this method we have run simulations 
with ^ = 0.6 eV where the number of TFs in the library were decreased from 500 by 
a factor of 5. On scales where the neutrinos contribute to the matter power spectrum, 
the difference found was below 0.05% for all redshifts. Also note that the non-linear 
neutrino correction term found is not due to an inaccurate sampling of the TFs, since 
such an error would manifest itself especially on the largest scales where the neutrino 
component contributes the most. 

We also tested our choices for the GADGET-2 parameters ErrTolIntAccuracy = 0.025 
and MaxSizeTimeStep = 0.03 by decreasing them by a factor of 4 and 6, respectively. 
At k < 1 hMpc~^ the difference is at most 0.1 % in each case, so regarding these two 
parameters the absolute power spectrum is evolved accurately in the region where we 
have the non-linear neutrino correction. We also ran the code with double precision and 
found an absolute error below 0.1%, which only manifested itself on small scales. 

Since the neutrino particle component mainly contributes to the gravitational 
potential in the region coverged by the PM grid, testing for convergence of the long 
range force is important, so as not to introduce errors when comparing the neutrino 
grid and particle methods. Therefore, we have run two simulations with a PM grid size 
of 1024 (models Cg and Cg). These simulations gave an almost identical calculation of 
the non-linear neutrino correction term compared to the A^pm = 512 case. 

In the standard TreePM approach the long range force is Gaussian filtered. Since 
we do not filter the neutrino grid modes, these modes contribute at smaller scales than 
do the neutrino particle long range modes. Since the PM approach breaks down at 
small scales, it is necessary to test how much this lack of filtering contributes. We have 
done this by running our 3 neutrino grid simulations (models Bi, Ci and Di) with the 
Gaussian filtering included. 

In the X^mi, = 1.2 eV case the non-linear neutrino correction term increases by 
only ~ 0.25% with the smoothing turned on. For the lower neutrino masses the effect is 
negligible (< 0.05%). Therefore, including the neutrino grid modes via the PM approach 
without a Gaussian filtering is well justified, since the neutrinos do not contribute 
substantially to the gravitational potential on scales where the PM method breaks down. 
Or equivalently, it is extremely well justified to neglect the linear neutrino component 
in the Tree part (for our chosen iV-body specific parameters, i.e. -Rbox = 512 h~^Mpc, 
NpM = 512, r, = 1.25 and r^ut = 4.5). 

4.4- The Extent in Fourier Space of the Neutrino Grid 

In order not to oversample the neutrinos in the grid approach it is usefull to find a 
maximum wavenumber beyond which the neutrinos do not contribute to the formation 
of structure. In the left panel of Fig. [H]we show the effect of neglecting the ^ m^, = 0.6 eV 
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Figure 8. Left: The effect on the square root of the matter power spectrum (an 
effective non-linear matter TF, Tm) of neglecting the perturbed neutrino component, 
= 0.6 eV, in the A^-body simulation when it is started at z = 49 (Cio/Ci). 
Right: The effect on Tm of neglecting the neutrino component, ^ m,^ = 0.6 eV, only 
when the power spectrum is calculated at a given redshift (model Ci). 



neutrino grid from z = 49 onwards, i.e. throughout the whole A^-body simulation. The 
neutrino component is still included in the calculation of the background evolution. The 
effect is of course pronounced at large scales, but even a lack of power is visible at small 
scales. The neutrino modes corresponding to these small scales do not contribute to this 
difference, as can be seen from the right panel of Fig. [HI Instead, it is power missing from 
a lack of mode-coupling between the small-scale CDM modes and large-scale CDM plus 
neutrino modes. Note that we have shown the difference in the square root of the matter 
power spectrum, since it is this quantity that is related to the gravitational potential 
calculated in Fourier space in gadget-2. 

To find a maximum wavenumber beyond which the neutrino grid does not 
contribute to the gravitational potential, the difference in the square root of the matter 
power spectrum with and without the neutrino grid is shown in the right panel of 
Fig. El for ^ rrii, = 0.6 eV. The grid is only taken out when the power spectrum is 
calculated and the neutrino modes are therefore included in the evolution of the matter 
perturbations. The wavenumber at which the difference goes to zero can then be used 
to find the required extent of the neutrino grid in Fourier space. For X^mi, = 0.6 eV 
a conservative value for this maximum wavenumber is ^ 1 /iMpc"^. Notice that this 
maximum wavenumber is not found by considering the power spectra today, but at a 
redshift of ~ 4. For '^m,^ = 0.3 eV the maximum wavenumber is ~ 0.5/iMpc~^. We 
caution that these maximum wavenumbers depend on Q^. 

5. Discussion and Conclusions 

We have presented a new method for implementing neutrinos in A^-body simulations 
which works extremely well for a total neutrino mass below ~ 0.5 eV. For such masses 
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the difference in the matter power spectrum compared with simulations where neutrinos 
are treated as particles is always below 1% on all scales. The precision is better for 
smaller masses, with the difference scaling roughly as ml, and the gain in computational 
speed compared to representing the neutrinos as A^-body particles is very large (scaling 
roughly as m~^). For all masses this is the only computationally feasible way to include 
neutrinos in simulations at the required level of precision, especially for high A^-body 
starting redshifts. 

We also note that the method presented here will work for any type of energy 
density which has almost linear perturbations. 

Finally we caution that although the method is very accurate for calculating all 
observables related to matter fluctuations, i.e. the power spectrum, halo mass functions 
etc, it is not accurate at the 1% level in describing the neutrino fluctuations alone. 
For predicting quantities such as the local density of relic neutrinos additional steps 
should be taken. One possibihty is to use the 1-particle Boltzmann technique [151 ITT] . 
In its simplest form, however, this method represents the neutrino component only in 
an approximate way. Another path is to solve with the method presented here until z 
becomes sufficiently low (in practise z < A) that the noise from the neutrino thermal 
velocity can be kept under control. At this point the grid simulation can be used as the 
initial condition for a simulation with neutrinos treated as particles. 
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